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Abstract 

In this work we study mixed mode oscillations in a model of secretion of GnRH (Gonadotropin Re- 
leasing Hormone) . The model is a phantom burster consisting of two feedforward coupled FitzHugh- 
Nagumo systems, with three time scales. The forcing system (Regulator) evolves on the slowest 
scale and acts by moving the slow null-cline of the forced system (Secretor). There are three modes 
of dynamics: pulsatility (transient relaxation oscillation), surge (quasi steady state) and small os- 
cillations related to the passage of the slow null-cline through a fold point of the fast null-cline. 
We derive a variety of reductions, taking advantage of the mentioned features of the system. We 
obtain two results; one on the local dynamics near the fold in the parameter regime corresponding 
to the presence of small oscillations and the other on the global dynamics, more specifically on 
the existence of an attracting limit cycle. Our local result is a rigorous characterization of small 
canards and sectors of rotation in the case of folded node with an additional time scale, a feature 
allowing for a clear geometric argument. The global result gives the existence of an attracting 
unique limit cycle, which, in some parameter regimes, remains attracting and unique even during 
passages through a canard explosion. 

Keywords : Slow-fast systems, multiple time scales, mixed mode oscillations, limit cycles, se- 
condary canards, sectors of rotations, folded node, singular perturbation, blow-up, GnRH secretion. 
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1 Introduction 



Mixed Mode Oscillations (MMOs) is a term used to describe trajectories that combine small oscil- 
lations and large oscillations of relaxation type, both recurring in an alternating manner. Recently 
there has been a lot of interest in MMOs that arise due to a generalized canard phenomenon, 
starting with the work of Milik, Szmolyan, Loeffelmann and Groeller [20]. Such MMOs arise in the 
context of slow- fast systems with at least two slow variables and with a folded critical manifold (set 
of equilibria of the fast system). The small oscillations arise during the passage of the trajectories 
near a fold, due to the presence of a so-called folded singularity. The dynamics near the folded 
singularity is transient, yet recurrent: the trajectories return to the neighborhood of the folded 
singularity by way of a global return mechanism. 

An important step on the way to an understanding of MMOs is the analysis of the flow near 
the folded singularities. Of particular importance are special solutions called canards. The term 
canard was first used to denote periodic solutions of the van der Pol equation that stayed close to 
the unstable slow manifold (approximated by the middle branch of the fast nullcline) [2]. One of 
the characteristic features of canard cycles is that they exist only for an exponentially small range 
of parameter values. This very sharp transition was then termed canard explosion [4]. The related 
term canard solution has been used to denote solutions connecting from a stable slow manifold 
to an unstable slow manifold. Such canards, sometimes also called maximal canards, organize 
the dynamics in a similar way as invariant sets which separate different dynamical regimes (e.g., 
separatrices of saddle points). In systems with more than one slow variable, canards occur in a 
more robust fashion and underlie the presence of the small oscillations near the folded singularity 
in MMOs. 

A prototypical example of a folded singularity with small oscillations is the folded node, studied 
by Benoit [1], by Wechselberger and Szmolyan [25], and by Wechselberger [27]. These articles 
focused on the local aspects of the dynamics. An exposition of how the dynamics near the folded 
node can be combined with a global return mechanism to lead to MMOs was given in [3] . This work 
was used as a basis of various explanations of MMO dynamics found in applications [24, 23, 11]. 
A shortcoming of the folded node approach is the lack of connection to a Hopf bifurcation, which 
seems to play a prominent role in many MMOs. This led to the interest in another, more degenerate 
folded singularity, known as Folded Saddle Node of type II (FSNII), originally introduced in [20] 
and recently analyzed in some detail by Krupa and Wechselberger [19]. Guckenheimer [12] studied 
a very similar problem in the parameter regime yet closer to the Hopf bifurcation, calling it singular 
Hopf bifurcation. For a more comprehensive overview we refer the reader to the recent review article 

[]■ 

Two notions that are central to the study of MMOs are secondary canards and sectors of 
rotation. Secondary canards [27, 3] are trajectories which originate in the attracting slow manifold, 
make a number of small oscillations in the fold region, and continue to the unstable slow manifold. 
There is ample numerical evidence of the existence and role of secondary canards ([8, 9]), as well 
as some partial theoretical results ([13, 27, 19, 16]). It has been highlighted that two trajectories 
crossing the region between two consecutive canards display the same number of small oscillations. 
Hence, the regions separated by secondary canards have been called sectors of rotation ([3]). As a 
parameter changes, a periodic orbit may move closer to a canard and pass to the adjacent sector 
of rotation. This transition has never been studied in detail. It is similar to a canard explosion, 
although more complicated, as chaotic behavior can be expected. 

The main result of this paper is that, in the context of our phantom burster problem, there 
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exists an attracting MMO orbit for all parameter values, also during the passage between different 
sectors. In addition we obtain a result on the existence of secondary canards of rotational type that 
is complementary to the results in [27, 16, 19] and relevant to the context of the phantom burster. 

It is important to note that, even if canards are more robust in three-dimensional slow-fast 
systems, they are still difficult to find numerically as well as particular types of MMOs. Forward 
integration is not possible due to the exponential expansion along the slow manifold, leading to 
exponential magnification of numerical errors. A breakthrough in the numerical detection and 
continuation of canards and MMOs has been achieved by Desroches et. al. [8] who used a boundary 
value approach in the context of numerical continuation with the software package Auto. 

In this article we investigate the presence of MMOs in the following system: 

eSx = -y + f(x), (la) 
5ij = a x + aiy + a 2 + cX, (lb) 

SX = -Y + g(X), (lc) 
Y = X + hY + b 2 , (Id) 

with 

/ : x — > A3X 3 + Aix, 
g : x -)• fi 3 x 3 + mx, 
A 3 ,/u 3 < 0, Xufii > 0, 
<H,c > 0, 
< e,S « 1. 

System (1) has been proposed in [5, 6, 26] to model the dynamics of GnRH secretion by hy- 
pothalamic neurons in female mammals. Subsystem (la)-(lb), called the Secretor, represents the 
mean-field approximation of the GnRH neuron population dynamics. It is driven, through the 
coupling term cX, by subsystem (lc)-(ld), called the Regulator, representing the activity of the 
interneuron population that conveys the periodic action of the ovarian steroids onto the GnRH 
neuron population. The time scale difference between the two oscillators is a transcription of the 
ratio between the ovarian cycle duration (few weeks) and the period of the secretory activity of the 
GnRH neuron population (few hours). 

From the point of view of dynamical classification, system (1) is a phantom burster with the 
additional feature that it has multiple time scales. Both the Regulator (lc)-(ld) and the driven 
Secretor (la)-(lb) are slow-fast systems of Fitzhugh-Nagumo type and there is a time scale difference 
between the two. The parameters controlling the time scales are e and S. 

We are interested in the case when the Regulator displays a stable relaxation limit cycle. Then, 
in a certain region of the parameter space (see [6]), the X-driven Secretor alternates between a fast 
oscillatory regime (when it displays an attracting relaxation limit cycle) and a stationary regime 
(when the current point tracks an attracting singular point). As a result, the signal generated by y, 
that represents the GnRH secretion along time, displays a periodic alternation of pulsatile regimes 
and surges as illustrated in Figure 1. During the transition from a surge back to the subsequent 
pulse phase, a pause corresponding to a segment of small oscillations may occur. 

In this article, we analyze the dynamical mechanism based on three different time scales that 
underlies the occurrence of the small oscillations. We prove that, for certain choices of the parameter 
values, the MMOs, including the pulse phase, surge and pause, exist and are stable limit cycles, 
even when close to a secondary canard. More precisely, we prove that canards with a specified 
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Figure 1: The signal y(t) generated along a typical orbit of system (1) displays a periodic alternation 
of pulsatile phases, surges and pauses. As displayed in the inset, the pauses consist of small 
oscillations. 

number of small oscillations are unique (with fixed choices of slow manifolds) and that any two 
adjacent canards differ by one rotation. Thus we prove that sectors of the same rotation (or simply 
sectors of rotation) exist and the passage, as a parameter varies, through a secondary canards adds 
(or subtracts) one small oscillation to the globally attracting orbit. 

The paper is organized as follows. In §2 we present the phantom burster dynamics of (1), discuss 
the different phases of the orbits and state the main result (Theorem 1). Section 3 is devoted to 
the local analysis of canard oscillations in a three time scale reduced system (Theorem 2). In §4 
we analyze the return mechanism and prove Theorem 1. Section 5 contains numerical findings 
illustrating our results. The article ends with a discussion section. 

2 Multiple time scale phantom bursting 

In this section we give a qualitative description of the dynamics we are interested in and state 
the main result. Parts of this section are a review and we refer the reader to [5, 6] for more 
details. We begin by sketching the basic features of the dynamics of the decoupled Regulator (lc)- 
(ld). Subsequently we set constraints on the Secretor's parameters in order to obtain the right 
dynamical behavior, introduce the different reduced systems suitable to describe various stages of 
the dynamics, and briefly describe the evolution of the x and y variables in the different stages of 
the dynamics. We end the section by stating the main theorem. 

2.1 The Regulator dynamics and its influence on the position of the Secretor 
slow nullcline and singular points 

We define 7 > by g' (±7) = so that the two knees of the cubic X-nullcline Y = g{X) are 
(±7, #(±7)). The knees split the cubic X-nullcline into three parts : the left and right branches 
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where g' < and the middle where g' > 0. As mentioned in the introduction, we assume that the 
parameters, specifically b\ and 62, are chosen so that the Regulator admits a relaxation limit cycle. 
To ensure this property, it is sufficient to assume that b\ is small enough (so that the Y-nullcline 
is steep enough) and that the y-nullcline intersects the cubic X-nullcline Y = g(X) on its middle 
branch (where g' < 0) away from the knees (see [(>]). Let us note X m \ a and X max respectively the 
minimal and maximal value of X along the Regulator limit cycle. For later reference we list four 
different phases of the evolution of X (see Figure 2): 

1. slow motion near the left branch of the cubic Y = g{X) : X increases slowly from X m [ Q to 

2. fast motion from the left knee to the right branch of the cubic : X increases quickly from —7 

to -Xmax> 

3. slow motion near the right branch of the cubic : X decreases slowly from X max to 7, 

4. fast motion from the right knee to the left branch of the cubic : X decreases quickly from 7 
to X m \ n . 




y = g(x) 

Figure 2: Limit cycle of (lc)-(ld) and its four different phases. 

The value of X drives the Secretor y-nullcline defined by a$x + a\y + 0,2 + cX = 0. Note 
that this nullcline is a straight line whose slope — ao/ai does not depend on X. As usual in the 
Fitzhugh-Nagumo system, a\ is assumed to be small so that the y-nullcline is very steep. 

As X increases (resp. decreases), the y-nullcline moves to the left (resp. to the right) in the 
Secretor phase space (x,y). Hence, depending on the value of X, the number of singular points 
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(lying on the cubic x-nullcline y = f(x)) varies. Also, their nature depends on their position with 
respect to the fold points (xf,f(xf)) and (—Xf,—f(xf)) that splits the x-nullcline into three part 
(left, middle and right branch). In particular: 

a. if a singular point lies on the middle branch outside a (9(e)-neighborhood of the folds, it is 
surrounded by a relaxation limit cycle ; 

b. if two different singular points lie on the left (resp. right) branch, the lowest (resp. highest) 
one is an attracting node and the highest (resp. lowest) one is a saddle ; 

c. if the Secretor admits a unique singular point, it is a saddle. 

The passage from a to b is a Hopf bifurcation that makes the limit cycle disappear through a canard 
explosion and the passage from b to c is a saddle-node bifurcation that makes the saddle and the 
node collapse. 

In the following, we assume (see hypotheses HI to H4 in the following section) that for all values 
of X between X m \ n and X max , the Secretor admits three different singular points determined by 
their x-component. Of special importance is the middle singular point (corresponding to the x-value 
lying between the two others) for which we note the x-component x s i ng (X). 

2.2 Constraints on the Secretor parameters and statement of the main result 

To obtain the qualitative behavior of the y-signal generated by full system (1) (Figure 1) we make 
the following hypotheses illustrated by Figure 3. 

(HI) The y-nullcline should pass through the right fold point of the cubic y = f(x) which generates 
the small oscillations. Hence, we assume that, for X = X mm , the y-nullcline should be on the 
right of - and close to - the upper fold (xf, f(xf): 

Xf X s [ n gyy\ m [ n J I.e. A m j n ^ v\f — 

This condition will be discussed in more detail later on. 

(H2) Once the y-nullcline has passed the right fold and the relaxation limit cycle of system (la)- 
(lb) appears, the cycle should persist until X = —7. Hence, we assume that for X = —7, the 
y-nullcline intersect the cubic y = f(x) on its middle branch: 

-x f < x sing (-7) i.e. - a x f + aif(-Xf) + a 2 - 07 < 0. 

(H3) From the beginning of the surge phase, system (la)-(lb) must admit an attracting node and 
a saddle on the left branch of the cubic y = f(x). This condition reads 

a + aiAi > 2v / aicAiX max 

which is equivalent to: 

Y ^ Y (ap + aiAi) 2 
4aicAi 

Let us note that value X$n of X corresponds to the saddle- node bifurcation of the Secretor 
occurring when the y-nullcline is tangent to the left branch of the cubic y = f(x). 



6 




7 



(H4) Until the end of the surge phase, system (la)-(lb) must admit an attracting node and a 
saddle on the left branch of the cubic y = f(x) as well. This condition reads: 

-x f > x sing (7) i.e. - a x f + aif(-xj) + a 2 + cj > 0. 

Figure 4 gives an instance of the relative positions of the Secretor slow nullcline that can arise 
due to the variation in X as (X, Y) traces the Regulator relaxation cycle under assumptions HI 
to H4. Then, the signal generated by variable y displays an alternation of surge, small oscillations 
and surge phases as illustrated in Figure 1. We refer to [ ] for an expanded explanation of the 
model behavior based on this approach. 




Figure 4: Instance of the Secretor slow nullcline locations under hypotheses (HI) to (H4). 
2.3 Reduced systems 

For each phase from 1 to 4 of the Regulator limit cycle, system (1) can be reduced using a specific 
approximation. We first recall the general process of desingularization applied to slow-fast systems 
near a fold. We introduce the different reduced systems that we will use for the dynamics analysis. 

Desingularized Reduced System. When a general slow-fast dynamical system 

ex = f(x,y,e), (2a) 
V = g(x,y,e), (2b) 
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is considered, with x and y of arbitrary dimension and < e -C 1, one classic way to understand 
the overall dynamics is by looking at the slow and the fast dynamics separately. An object of 
great importance for both the slow and the fast dynamics of the full system is the so-called critical 
manifold C° defined as the nullcline for the fast variable, that is: 

C° = {(x,y); /(s,y,0) = o}. 

Consequently, the critical manifold is the phase space of the reduced system obtained by setting 
e = in equations (2) and which approximates the slow dynamics of the original system; the 
reduced system is a differential-algebraic equation. In order to understand the flow of the reduced 
system, which then takes place on C° and is associated with the singular limit e = 0, the usual 
strategy - which we will use several times in the rest of the paper - is to differentiate the algebraic 
equation defining C° with respect to time. This gives 

xf x (x,y,0) + yf y (x,y,0) = 0, 

V = g(x,y,0), 

which reduces to 

xf x (x,y,0) = -g(x,y,0)f y (x,y,0), (3a) 

y = g(x,y,0). (3b) 

The previous system is singular along the fold set of C° with respect to the fast variable x, that 
is, the set T = {f x (x,y,0) = 0}. In order to understand the slow flow up to the fold set, one can 
desingularize system (3) via a rescaling of factor f x (x, y, 0), which yields the so-called desingularized 
reduced system. In that case, special care has to be taken going from the desingularized reduced 
system back to the slow system. Indeed the previous rescaling changes the orientation of orbits 
when f x (x,y,0) < and one needs to reverse orientation in order to get the correct direction of 
the flow in the original reduced system. 



Three dimensional reduction with three time scales during the pulsatile phase. Under 
the preceding assumptions, slow motion 1 (X m ; n < X < 7) corresponds for system (la)-(lb) to the 
oscillatory phase producing the small oscillations and subsequently the pulses in the y-signal. 

The variables X and Y follow the slowest time scale and the current point (X, Y) remains in 
a 0(<5)-neighborhood of the cubic. This reads Y = hs(X) where (X,S) \-t h$(X) is an analytic 
function on ] — 00, — 7[xM^_ and h = g. Thus, on ] — 00, — 7[, h' s (X) = g'(X) + 0(5). 

We introduce a reduced system obtained from (1) assuming that Y = hs(X). We differentiate 
this condition, with 5 constant: 

Y = Xh' 5 (X). 

By replacing the dynamics of Y in (1), one obtains the three-dimensional system with three different 
time scales: 

-y + f(x), (4a) 
a x + a\y + a 2 + cX, (4b) 
X + b 1 (g(X) + 0(5)) + b 2 

g'(X) + 0(5) ■ [ ) 



e5x = 
Sy = 

X = 



9 



Two-dimensional reduction with two time scales during the surge phase. Slow motion 
3 (7 < X < X max ) corresponds to the surge phase. The current point (x,y) follows the attracting 
node of (la)-(lb) lying on the left branch of y = f(x). Hence, both approximation X ~ g(X) and 
y ~ f(x) stand. 

By reducing the fastest time scale, i.e., by setting y = f(x) in (1), we obtain the following 
system: 

5f'(x)x = aox + a\f(x) + a 2 + cX, (5a) 
SX = -Y + g(X), (5b) 
Y = X + hY + b 2 . (5c) 

Then setting Y = g(X) leads to the following equations: 

Sf (x)x = aox + a\f(x) + a 2 + cX, (6a) 
g'(X)X = X + b l9 {X) + b 2 . (6b) 

Away from the folds of both cubics (where f'(x) = or g'(X) = 0) we can rewrite (6) as follows: 

x . _ a x + aif(x) + a 2 + cX 

dX ~ f{x) ' (7&) 

■ X + b ig (X) + b2 , . 

x - gvn • (7b) 

Hence we have obtained a two-dimensional slow- fast system with slow variable X and fast variable x. 

Boundary-layer system during the transitions. During fast motions 2 and 4, {X, Y) evolves 

according to the X time scale and the slowest variable Y is almost constant. By setting 5 = in 
the rescaled system 

ex = -y + f(x), (8a) 

y = a x + aiy + a 2 + cX, (8b) 

X = -Y + g(X), (8c) 

Y = 5(X + b x Y + 62), (8d) 

one obtains the so-called Boundary-Layer System 

ex = -y + f(x), (9a) 
y = aox + aiy + a 2 + cX, (9b) 
X = -Y + g(X), (9c) 

where Y ~ g(— 7) for fast motion 2 and Y — 9(7) for fast motion 4. 
2.4 Different dynamical regimes 

Surge. The surge corresponds to phase 3 of §2.1, when X passes from X max to 7. The dynamics is 
governed by system (7). Initially x decreases to reach the vicinity of the nullcline X = f(x), where 

~ a x + a 1 f(x) + a 2 
fix) = ; • 
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This is coupled with a significant increase of y. Subsequently x grows and y decreases, moving at 
the rate given by the slowest time scale 0(1/ (5e)). 

Small oscillations during the post-surge pause. Hypothesis (HI) guarantees that the surge is 
followed by a sequence of small oscillations taking place near the fold (xt, f(xf)). As X has already 
reached the vicinity of X m \ n the dynamics is governed by (4) and can be described as follows. After 
the surge, the trajectory is attracted to a stable quasi steady state of node or focus type on the 
right branch of y = f(x). As X increases from A" m ; n , the quasi steady state changes stability as 
a pair of complex eigenvalues passes through the imaginary axis, that is, a slow passage through 
a Hopf bifurcation occurs [21, 22]. This phenomenon constitutes a delayed transition from surge 
to pulsatility. The analysis of the small oscillations (see Figure 5) constitutes a large part of this 
article. 




Figure 5: Small oscillations corresponding to the post-surge pause. Panels (a) and (b) the full 
periodic orbit and a zoom in the region of the pause, respectively, in the three-dimensional phase 
space (x,y,X); panel (c) shows a stronger zoom on the small oscillations of the pause, together 
with the critical manifold C°, the fold curve J- (red line) and the folded node (dot) of the three- 
dimensional subsystem (4). 
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Recall that Xf (introduced in hypothesis (HI)) is defined as the value of X for which the y- 
nullcline intersects the cubic a>nullcline at its right fold. We will assume that X m i n < Xf for the 
following reason. If X m [ Q is significantly less than Xf then the quasi steady state is still a node 
at the moment when the trajectory is attracted to it. During the passage near the fold the quasi 
steady-state turns into a stable focus and subsequently becomes unstable; first it is an unstable 
focus and then an unstable node. However, the trajectory is extremely close to the quasi steady 
state when it is a focus and only gets repelled from it after it has changed from unstable focus to 
unstable node. Hence no small oscillations can be seen. On the other hand, if X m ; n < Xf then the 
aforementioned quasi steady state can always be a focus, initially stable and subsequently unstable. 

An important aspect of the dynamics are canards. A canard segment is a segment of a trajec- 
tory which initially stays close to a stable branch of y = f(x) for a time of 0(1 /e), subsequently 
passes through the fold region and finally remains near the middle branch of y = f(x) for a time 
of 0(l/e). A canard is a trajectory containing a canard segment. When a system possesses a 
folded singularity, there can be canard trajectories with small oscillations. We then define a k-th 
secondary canard as the canard trajectory making k small oscillations near the fold. As part of the 
analysis we show that secondary canards separate the trajectories with different numbers of small 
oscillations. 

Pulsatility. Pulsatility is a region of transient relaxation oscillation corresponding to phase 1 of 
§2.1. It is a direct continuation of the pause and the governing system is still (4). As the slow null- 
cline of (la)-(lb) cuts through the middle branch of y = f(x) the dynamics is purely of relaxation 
type. 

Transitions. The first transition corresponds to phase 4 of §2.1 and follows the surge. The 
governing system is (9) with Y = 5(7). The variables (x,y) first evolve on the intermediate 
time scale 0(1/ e), following the nullcline y = f(x), subsequently jump to the right branch of the 
nullcline y = f(x) and then follow the right branch of y = f(x), evolving on the time scale 0(l/e) 
and arriving to the vicinity of (xf,yf) as X reaches the vicinity of X m \ n . 

The second transition corresponds to phase 2 of §2.1 and precedes the surge. The governing 
system is (9) with Y = g(— 7). Before X reaches the vicinity of AT max , the point (x,y) can move 
down along the left branch of y = f(x) and then turn back up the left branch of y = f(x) or 
it can jump to the right branch of y = f(x) making another pulse before the surge. There is a 
canard phenomenon associated with this behavior which can incur some expansion. Estimating 
this expansion is a part of our analysis. 

2.5 Statement of the main theorem 

Before we can state our main theorem, we need to make one additional assumption: 




2 

> - 

e 




where 



/(*) 



aox + aif(x) + a 2 



c 



Hypothesis (H5) guarantees that the return map around the cycle is contracting. 
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Theorm 1 Provided that 5 and e are sufficiently small and (H1)-(H5) hold, there exists a unique 
stable limit cycle consisting of a number of small oscillations, a number of pulses and one surge. 
Some exceptional limit cycles, existing only in exponentially small parameter regions, contain canard 
segments. All the limit cycles are fixed points of a single passage around the cycle of surge, pause 
and pulsatility. Varying a regular parameter can lead to a change in the number of pulses or small 
oscillations by means of a passage through a canard explosion. There are two canard explosions, one 
associated with the upper fold and one with the lower fold. A passage through the canard explosion 
at the upper fold yields a transformation of a small oscillation to a pulse or vice versa. The passage 
through the canard explosion at the lower fold leads to an addition or a subtraction of a pulse. 



3 Folded singularities of system (4) 

Folded singularities are usually studied in systems with two slow variables, however in this paper 
we need to consider system (4), which has three time scales. The usual approach for classifying 
folded singularities is to consider the desingularized reduced system, see [25] for instance. Here we 
will mimic this procedure for (4). 



3.1 Nature of the folded singularity 

Since the fastest variable in (4) is x, we set the left hand side of (4a) to 0, obtaining the constraint 
y = f(x). Hence the critical manifold corresponding to the fastest time scale is the cubic surface 
C° = {y = f{x)}. It displays two folds respectively for 



±x f = ± ] J^-, y = f{x). 

By applying the procedure described at the beginning of §2.3, one obtains the desingularized reduced 
system 

x = — {aox + a\f(x)+a 2 + cX), (10a) 
• / X + b 1 (g(X) + 0(8)) + b 2 \ 

x = gl{x) - 0{§) jf(x) = e 5 (x,x). (iob) 

Note that the slow flow (4b)- (4c) on C° has the same orbits as the desingularized reduced system 
(10), however one has to reverse the orientation where f'(x) > 0, that is, on the repelling sheet of 
the critical manifold C°. 

The equilibria of system (10) on the fold curve J 7 , that is, the folded singularities of the 3D 
system, are given by 



f'(x f ) = 0, (11a 
X f = --(a Xf + aiy f + a 2 ) 
with yf given by yf = f{xf). From the expression of f(x), we get 



Xf = -^(a x f + aiy f + a 2 ), (lib) 



*/ = (12a) 
X f = -lU^-(a + ^X l a 1 ) + a 2 ), (12b) 



13 



with 



2 / -A 

Vf = x f( X 3x} + Ai) = igAi^ 



Note that 

due to the factor f'(x). In addition we have f"(xf) = 6A3XJ. Hence, the jacobian matrix j( 10 ) of 
system (10) at (xf,Xf) reads 

•*"»<*'•*'> = (-easx, "oj- < 13 > 

The eigenvalues of the matrix J^ 10 ) are given by 



It follows that if 

X eval = 2AcX 3 x f ( Xf + b g ) 9 { ^ + b2 ) < (15) 

then, for small enough S, there are two real eigenvalues of the same sign, i.e. the folded singularity 
is a folded node. Evaluating (15) for c = 0.69, ao = 1, a\ = 0.02, a 2 = 0.8, b\ = 0, b 2 = —0.8, 
A3 = —1, Ai = 1.5, fi3 = —1, m = 4, one obtains X eva i = —3.3248. 

3.2 Local form near the folded singularity 

We translate the origin to (x f, f(xj), X f), with (xf,Xf) given by (11), and rescale the x and y 
variables. We first set 

x = a(x — xj) , 
V = a>{y-f(x f )), 

where a = \/— 3A1A3. In these new coordinates, system (4) reads 

ebx = -y-x - —x , 

5y = a x + aiy + a(a Xf + aif(xf) + a 2 + cX) , 
> 1 (g(X) + Q(5) 
g'{X) + 0{5) 

Now we translate the variable X by setting X = X — Xf to obtain the following system 

s5x = — y — x 2 — — — x 3 , (16a) 
3Ai 

Sy = aox + a%y + acX, (16b) 

X = ip + TpX + 0(X 2 ,5). (16c) 



^ X + b 1 (g(X) + 0(S)) + b 2 
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where 

X f + hg(X f ) + b 2 
V = Ti — \ • 

Finally we introduce new variables (x,y,z) denned by (z, — x,y) = (x,y,X) and rescale the time 
to obtain the system: 

x = acy — aoz + O(x), (17a) 
y = 5{^ + ^y + 0(y 2 ,5)) , (17b) 
ez = x + z 2 + 0(z 3 ). (17c) 

For 5=1, system (17) is analogous to the normal forms considered in [J] and [25]. Assuming that 
e and 5 are small constants (perturbation parameters), this system has one fast, one slow and one 
super-slow variable. Note that X eva i = 24cXsXjp, A3 < and c,Xf > 0. Hence, if (15) holds then 
(p > and the folded singularity is a folded node. The case if < 0, or equivalently X eya \ > 0, gives 
a folded saddle for 5 sufficiently small. Finally <p = corresponds to a folded saddle-node. Other 
types of folded singularities do not occur for 8 close to 0. 



3.3 Local analysis near the folded node: statement of the result 

After dropping the ~ signs and rescaling time, system (16) reads: 

n 1 q , . 

ex = — y — x — ttt-x , (18a) 
3Ai 

y = aox + a±y + acX, (18b) 

X = 5{ip + ^X + 0(X 2 ,5)). (18c) 

Let the section S m defined by y = —p 2 , where p > is small but fixed. Let S a , £ be the 
attracting Fenichel slow manifold, perturbed from the critical manifold y = x 2 + x 3 /(3Ai), near 
the section S m . Note that E m is a transverse section of the flow of (18) intersecting S a e close to 
the fold but O(l) away from it. In this section we focus on describing the dynamics starting in the 
curve S at£ n S m . Each trajectory starting in S a<£ n S in enters the neighborhood of the fold, makes 
a number of small rotations, and then exits the fold region. The number of small oscillations can 
be different for different trajectories. Canards can now be defined as the trajectories that go into 
the repelling slow manifold SV j£ (see §2.4 for an alternative definition). A k th secondary canard is 
a canard that makes k small oscillations in the fold region and subsequently runs into S rjE . 

Suppose the number of the small rotations for two trajectories (x, y, X) and (x, y, X) is different. 
Then there exists a secondary canard with initial condition somewhere on the segment between 
(x,y,X) and (x,y,Y). This way we can define sectors of the same rotation, or simply sectors of 
rotation, as the segments of S aiE n S m between the consecutive canards. We now state the main 
theorem of this section. This theorem leads to a precise definition and description of the sectors of 
rotation. 



Theorm 2 There exists a number R > such that, for every < v < R there exists a family of 
k th secondary canards with 

v , R 
-5 <k< J- 
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The canards with consecutive rotation numbers are next to each other. The distance between the 
consecutive canards measured in the section E in is bounded below by G\b~\fe and above by Cib^fe, 
where C\ and C2 are positive constants. 

Corollary 1 The k th sector of rotation, defined as the region between the k th and the (k + l) s * 
secondary canard consists of points whose trajectories make k rotations in the fold region. 

Our proof of Theorem 2 is based on the application of the blow-up method and builds on the results 
of [19]. This section is organized as follows. In 3.4 we introduce the blow-up and its charts. In 3.5 
we analyze the dynamics in the entry chart Kl, which covers a region near the critical manifold 
from S in to 0(y/e) away from the fold. In 3.6 we analyze the dynamics in the central chart K2, 
which describes the region very close to the singularity. We describe the delayed Hopf bifurcation 
occurring in this chart relying on the results of [ 9]. In 3.7 we prove Theorem 2, building on the 
results of Sections 3.5 and 3.6. 



3.4 Blow-up 

We use the following blow-up function: 



$ : l + x5 4 ^ M 5 , 

(f, x, y, X, e) -> (fx, f 2 y, fX, f 2 e) = (x, y, X, e). 



(19) 



In the entry chart y = —1, the blow-up (19) transforms variables (x,y,X,e) into new variables 
(x\, ri, Xi, £1) as follows: 

x = rixi, y = -r\, X = r 1 X 1 , e = r\e\. (20) 

After the transformation of system (18) and omitting a factor of r\ (which corresponds to a time 
rescaling), we obtain the following equations: 

xi = --xie x F(xi,ri,Xi) - (-1 + x 2 +rx—-x\), (21a) 

I OAl 

fi = -riEiF(xi,ri,X x ), (21b) 

Xi = -^Xi6iF(xi,n,si) + 6iS(cp + 7PnXi + 0(rfxlS)), (21c) 
ii = -e\F(xi,ri,Xi), (21d) 

where F(xi,ri,X\) = —a$xi + a\r\ — acX\. 

In the transition chart e = 1, the blow-up (19) corresponds to the change of variables (x, y, X, (p, e) 
into {x2,y2,X 2 ,(p2,r2) defined by 

x = r 2 x 2 , y = rjy 2 , X = r 2 X 2 , e = r\, (22) 

that is equivalent to 

x = *fex 2 , y = ey 2 , X = ^feX 2 . (23) 
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After the transformation of system (18), canceling a factor of yfe and canceling the equation e = 0, 
we obtain the following equations: 



x 2 = -2/2 - x% - ( 24a ) 
oAi 

2/2 = 00x2 + 00X2 + 01^2/2, (24b) 

X 2 = 5{^ + ^X 2 + 0{eXl,5)) . (24c) 

To prove Theorem 2, we need to find trajectories connecting from S a ,e to S r , e - As it will become 
clear from our forthcoming analysis, there is a natural extension of S rE to Kl. Hence we will be 
following the dynamics from Kl to K2 and then back to Kl. Consequently, we need to be able to 
transform Kl to K2 and vice versa on the overlap of the charts. The transformations between Kl 
and K2 are given by 

x 2 = ^, 2/2 = --, X 2 = ^=, (KltoK2), (25) 
v^i £1 V £ l 



and 

xi = -^=, ei = -— , Xi = -^=, (K2to^l). (26) 
V-2/2 2/2 V-2/2 

3.5 Extending Fenichel theory in chart fQ 

The key observation concerning the dynamics of (21) is that there exist center manifolds defined, 
approximately, by x\ ~ ±1. These manifolds are the extensions of the Fenichel slow manifolds S a}£ 
and S r>e . Near r\ = 0, they intersect the hyperplane E\ = according to the equation 

-1 + xi 2 + n^-x\ = 0. 

The center manifold CM a corresponding to x\ ~ 1 is attracting and given, near ri = 0, = 0, by 
the development 

xi = 1 - ^r-ri + 0(si,rf). 

Similarly, there exists an unstable center manifold CM r , corresponding to x\ ~ —1. 

The restriction of the flow to CM a , after canceling a factor of e%, which amounts to a time 
rescaling, is given by 

h = -riF(ri, (27a) 

Xt = -^X 1 F{r 1 ,X 1 ,e 1 ) + 6{ip + iPr 1 X 1 + 0{r 2 l Xl6)), (27b) 
ex = -eiF(n,Xi,ei), (27c) 

with 

F(n, Xi,ei) = -00 + (ai + ^r-)n - acA"i + 0(ei,r?) 

DAi 

being the restriction of F to CM a . 
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Hyperplanes r\ = and E\ = are invariant for (27). As if > 0, this system admits two singular 
points lying in r\ = s\ = and denned by their X\ component: 



± -^/oq — 8<5ac(£> 



2ac 



(28) 



Both values of X\ are negative and, since 5 is small 

a 



ac 



< x- < xt < o. 



Both singular points are hyperbolic. 

At each of the singular point (ri, X\, £i) = (0, X x , 0), the jacobian matrix J\ associated with 
system (27) reads 



jl 27) (0,Af,0) 



V 



-\{a + acXf) 

■lXt(ai + ^j+5{i>Xt + 0(S)) \{a + 2acXf) -\<D(l) 
a + acX^ 



For 5 — > + , one obtains at (0, X^ 1 ", 0), 



J ( Sf ) (0,X+,0)= I 0+ a -\o(i) I , 



2 u 







and, at (0, X l , 0), 



j ( f (o,xr 



0) 



V 



2ac 



«1 + 




go 
6Ai 



2 








For 5 > small there exists a two-dimensional center manifold associated with the equilibrium 
(0,X[~,0). For 5 = this manifold is defined by the condition F(r\, Xi,e\) = and consists 
entirely of equilibria. For 5 > the flow on the center manifold becomes weakly hyperbolic. To 
estimate the corresponding eigenvalues, we use (28) to obtain the approximation 



X r = -^ + ^o + 0( 2 ). 

ac ao 



(29) 



Using the fact that ao, a, c and (p are positive it is easy to see that the eigenvalue corresponding to 
the ri-direction is negative and the eigenvalue corresponding to the ei-direction is positive. Hence 
r\ decreases and S\ increases along trajectories. The flow of (27) is shown in Figure 6 and will be 
discussed in more detail below. 

Based on the information collected above we can follow the passage of CM a , and thus S aiE , 
until the entry in the fold region. We define the following sections of the flow of (27a)-(27c): 



{(x 1 ,r 1 ,X 1 ,e 1 ,ip 1 ) : n = p}, 
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Figure 6: Flow of (27) obtained by reduction of the flow of the blown-up system in chart Kl to the 
attracting center manifold CM a . It induces a transition between the sections S 1 " and S° ut illustrate 
by the red orbits. The repulsive two-dimensional center manifold (sky blue) associated with the 
singular point (0, X^, 0) (near the origin) contains the one-dimensional stable manifold W S {X^) of 
this singular point the strong unstable manifold W uu (Xf ). The attracting center manifold (purple) 
associated with the singular point (0,X^,0) is built in the same manner. Along the flow, the X\ 
direction is expanded away from the sky blue manifold and contracted towards the purple one. 

where p > is the constant used in the definition of E m , and 

T% at = {(x 1 ,r 1 ,X 1 ,e 1 ,ip 1 ) : e x = t?}, 

where rj > is a sufficiently small constant. We begin by explaining the meaning of E 1 " and £° ut in 
the context of the coordinates (x,y,X) of system (18) and (x2,U2, X2) of system (24). The section 
E^ n corresponds to the section S m defined at the beginning of §3.3. Further it follows from (25) 
that the section E° ut transforms to the section 

4 n = l(x 2 ,y2,X 2 ) : V2 = ~^\- 

Recall that S atE is the attracting Fenichel slow manifold near the section S m . The set S a ^ £ n S m is 
approximated by the line 

{(x,y,X) : x = p,y = -p 2 } . 
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To see which points in S° ut can be reached by trajectories starting in CM a n Xf, we study the 
dynamics of system (27) whose flow approximates the flow of (21) on CM a . We see that the X±- 
direction is expanded away from the equilibrium at the origin and is contracted towards the center 
manifold of the equilibrium (0,X^ ,0). Hence the projection of the center manifold CM a onto the 
Xi-direction contains the interval (— ao/(ac), 0). In fact the intersection of the center manifold 
with S° ut is a thin band containing a line segment which is close to the interval 

£i = l, -— <Xi<0. 

ac 

3.6 Delayed Hopf bifurcation and the way in/way out function in chart K2 

We use system (24) to compute the way-in/ way-out function near the fold. Note that (24) is a 
slow-fast system with two fast and one slow variables, with singular parameter 5. The critical 
manifold of (24) is given by 



\ 2 



2 



ac ac 

X2 = A 2 , J/2 = - — A 2 • 

a \a J 

We will denote this manifold by So- The linearization of the fast system about So is given by the 
matrix 

2A(e)X -1 



(I ' (30) 
ao U 

with A(e) = ac/ao + 0(\/e). Note that A(e) > 0, for sufficiently small e. Hence, by Fenichel 
theory, given a constant 77 > 0, there exist slow manifolds S- g and S +i s, attracting and repelling 
respectively, close to the line segments 

-- < X < -7] and 77 < X < - 
77 7/ 

respectively. 

We visualize the flow in the original coordinates (X2, y2,X 2 ) with X2 < 0, in Figure 7 ; the flow 
is strongly contracting in a tube around the critical manifold So ■ Similarly, Figure 9 shows the flow 
near So with X2 > 0; there the flow is strongly expanding. 

We rectify So by translating it to the line (0,0,^2), which is achieved by a transformation of 
the form: ^ 

x 2 = x-— X 2 + 0(v^), y 2 = y-( aC ) Xf + O^), X 2 = X. 
ao \aoJ 

In the new variables, system (24) reads (after dropping the 7 signs) 

x = -y + 2A(e)Xx- x 2 + 0(y/e,6), (31a) 
y = a x + 0(y/e, 5,x 2 ), (31b) 
X = 5((p + ^X + 0( v / iX 2 ))., (31c) 

Note that y/e is a regular parameter in (31). Hence, to simplify the notation, we will suppress the 
dependance of A on e. In Figure 8 we show the flow in the rectified coordinates for both X negative 
and positive. 
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f 2/2 




Figure 7: Transition from Erf to the cylinder Cs in chart -fT2. The image of 5 a)E nS m in E™ (purple 
interval) lies in X2 = 1/ y/rj. Under the flow of system (24), trajectories starting from the purple 
interval reach the vicinity of the critical manifold Sq by entering the cylinder C$ . 



Note that the eigenvalues of the matrix given by (30) are off the real axis if 



a 0<x< V«o 



A A 

with negative real part for X < and with positive real part for X > 0. For X < 0, we define the 
function ^f(X) by the formula 

*(X) = X* with I — - — dZ = 0. (32) 
Jx <P + ^Z 

The function ^ is the way in/way out function for all X satisfying 



an 

< X < 0. 



A 

Heuristically this means that the trajectories attracted to S-^ near X, will be repelled from 
near X*. To state a more precise result we introduce, for Xq < 0, the sections 

S Xo ={(x,y,X) : X = X }. 

The following result characterizes the transition map from a section Sx to Sx* ; ^* = ^(Xq). 
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Proposition 1 There exist constants r\\ > and rj2 > with the following property. For any 



«0 



A 



and (x,y,Xo) £ T,x sufficiently close to the origin in T,x , let (x(t),y(t), X(t)) be the trajectory 
of (31) starting at (x,y,Xo). Let X* = ^(Xq). There exists t* > such that 

X(t*) = X*, 

Vt G]0,t*[, X < X(t) < X* and x(t) 2 + y(t) 2 < 771. 

Moreover, if the distance between (x,y,Xo) and S-^ n T,x equals £ 2 ( 1+Q! ) /or some < a < 1/4 
i/ien t/ie distance between X*) and S4. 5 is bounded below by 772<5 2 ^ 1+c ^. 

A proof of this result based on the work of Neishtadt [22] can be found in [19], Corollary 5.1. 
Reference [19] also contains extensions of Proposition 1 to the case of Xq < y/a~Q/A, but we will 
not be concerned with these results in this article. 



3.7 Proof of Theorem 2 

In this section we will use both the original coordinates (X2, 2/2,-^2) ano - * ne rectified coordinates 
(x,y,X) in the following way. Let C$ denote a cylinder of radius <5 2 ( 1+Q ) around S- t s, where 
< a < 1/4 is a constant. Trajectories starting in S 2 n , with X2 < will enter C5 (for small 
enough 5) and subsequently exit C§, with X2 > 0. We will use the (x%, y 2 ,X 2 ) for the part of the 
trajectories outside Cs and (x,y,X) for the part of the trajectories inside C$. 
Recall that the interval 

xi = 1, - — < Xi < 
ac 

is the image of S a ,e H S m in E° . Translated to S 2 n this interval has the form 

*2 = ^= < X 2 < 0. 

Recall the center manifold CM a which coincided with the extension to K\ of the slow manifold 
S a ,e- We denote the image of CM a in K2 by transformation (25) also by CM a . Note that the point 



{x 2 ,y 2 ,X 2 ) 



1 1 a 



y/fj rj ac^j 



is on the critical manifold of system (24). By choosing r] < 1/oq, we guarantee that the point where 
the eigenvalues change from real to complex, which is given by 



X2 = 

is included in the interval (— ao / (ctCy/rj) , 0) provided that 5 is small enough. Consider the segment 
of trajectory of (24) starting at a point in CM a n S 2 n suc h that 

< X < -r/o < 



A 
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and ending at a point in Cg. Note that the flow is predominantly in the fast directions and away 
from singularities. Hence the passage time is 0(ln<5) uniformly in 5 and in X2 6 (—^/a /A, rjo). It 
follows that the ^-coordinate of the endpoint of the segment of trajectory also satisfies 

-^<x 2 < m , 

provided that 5 is small enough. 

We proceed using the rectified coordinates (x,y,X). Consider a trajectory starting at a point 
in C$ with X = Xq. Such trajectory first follows S-j and subsequently S +! $ until X becomes 
approximately equal to ^f(X), see (32). Corollary 5.1 in [ 9] states that the flow of (24) for 
trajectories starting at Cg is linear at lowest order (this occurs for the specified choice of a). 
Moreover, for each Xq verifying 

V "° <X < -770 < 0, 



A 

there exists X*(Xq, 5) such that the transition from Cg fl {X = Xq} to {X = X*(Xq, 5)} with 

X*(X ,5) =X*(X ) + O(5) 
is, at lowest order, a pure rotation by the angle R(Xq)/5, where 



Jx f + tpZ 

We refer the reader to [19], §5, for further detail. 

Consider an interval (xq, yo, X) G Cg with X varying within 0(5) of Xq and for every X consider 
the image of (xq, yo, X) by the transition via the flow of (24) to the section {X = X*(Xq, 5)}. The 
image of the interval is a segment of a very tight spiral, as changing X by 0(5) produces an 
increment of the angle greater than 2tt while X*(Xo,6) changes very little. Now consider the 
continuation of CM r backwards in time, near X*(Xq, 5). Since, backwards in time, the trajectories 
on CM r follow closely the fast fibers of (24) and converge to £+ 5, they must intersect the mentioned 
spiral transversely. If, instead of taking the interval {(xo,yo, X) G Cg,X = Xq + 0(5)}, we take 
a segment of the continuation of CM a to Cg, we obtain a similar spiral and similar transverse 
intersections, separated by a distance bounded below by K5, for some constant K > (see Figure 
8). 

These intersections correspond to secondary canards. A computation shows that 



A 



< X <-r ] ^ V'(X ) < 0. 



Hence the number of rotations of the secondary canards monotonically increases as Xq decreases. 
This implies that secondary canards are unique and that the consecutive canards differ by one 
rotation. It also follows that the distance between the secondary canards is 0(5). After translating 
back to the original coordinates (x,y,X) (blowing down), we obtain the estimate of Theorem 2. 
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Figure 8: Way-in/way-out transition near the rectified critical manifold of (31). The trajectory 
entering the cylinder Cs at X = Xq < is first attracted by the slow manifold S- t $, remain near 
5 +j 5 for a while and exits the cylinder at X = ^(Xq) > 0. A 2ir increment of the exit angle in 
(x, y) is obtained by a 0(<5) variation of the value Xq of the entry. 




Figure 9: Transition from the cylinder Cg in chart K2 to S re H XJf (purple interval). The figure 
shows segments of trajectories of (24) starting in the interior of the cylinder Cs and ending in 
S r>£ n £!f. This figure can be understood as the 'backwards in time' version of Figure 7. 
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3.8 Passage through the fold region — composition of the dynamics in the charts 

Using Theorem 2 we can now describe the dynamics in the fold region, starting from the union 
of the sectors of rotation in the section S in to the exit from the fold region. We will restrict our 
attention to the trajectories on S a<e as all the other trajectories shadow a trajectory on S a e . We 
first consider the trajectories that are not close to canards. After passing through C$ and exiting 
through its boundary for X > these trajectories separate quickly from 5 r>e and exit along the fast 
fibers. This is a simple and fast transition which does not incur much contraction. This transition 
is difficult to study mathematically due to resonance, but we will not focus on it here, referring the 
reader to [17]. 

If the trajectories are close to a canard they will arrive in S m and must continue to either E* , 
which is defined by x 2 = 0, (in the original coordinates by x = Xf), and subsequently reach S in , or 
they pass to the left branch of the nullcline resembling a canard with head and subsequently reach 



4 Proof of Theorem 1 

4.1 Contraction during the surge 

Recall that surge begins as X approaches X max and the dynamics is governed by (7). We rewrite 
(7) for convenience using a different notation: 

X - f{x) 

= ° f'(x) ' ( ' 

■ _ x + Mi) + b 2 

x ~ Tlx) ' (33b) 

with 

~. . a x + aif(x) + o 2 -aiA 3 x 3 - (a + aiAi)x - a 2 

f{x) = = ■ 

c c 

Recall the definition of x S i ng following the statement of Theorem 1 and note that 



max • 



/(3'sing(A max )) — X 
/0sing(7)) = 7- 

The slow manifold of (33) is defined by X = f(x). It turns out that / is an S shaped curve. To 
show that we compute the critical points. There are two, given by the formula: 

x 2 = - a ° + a x lAl ^ 30. (34) 
3aiA 3 

Further, f"(x) = — 60^32;, hence the negative critical point is a maximum and the positive one a 
minimum. Let x c + > > x c - be the critical points. Note that hypothesis (H3) is equivalent to the 
condition X max < f(x c J). The slow and fast dynamics of (33) are shown in Figure 10. 

It follows from the slow-fast structure of (33) and from the hypotheses (H3) and (H4) that 
the minimal (resp. maximal) value of x during the surge is close to x s i ng (X max ) (resp. x smg (7)). 
The passage through surge is always an exponentially strong contraction, with contraction rate 
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Figure 10: Slow and fast dynamics of (33) that approximate the flow of the whole system during 
the surge. Hypothesis (H3) ensure that X is always smaller than the value of X at the knee. Hence, 
at the beginning of the surge, x is near and X = X max . Under the fast dynamics (33a), (x,X) 
quickly reaches the slow manifold vicinity and then slowly goes down along it following the slow 
dynamics (33b). This mechanism is known to induce an exponential contraction between orbits. 



0(exp(— C3/5)), where C3 > is a constant. To compute C3 we follow the approach of [18], §2. 
We consider the reduced problem of (33), given by: 





X 



or, parametrized by x, 



= -f(x)+X, 

X + b l9 (X) + b 2 
9'(X) 

f(x) + b ig (f(x)) + b2 



(35a) 
(35b) 

(36) 



g'(f(x))f'(x) 

Let xo(t) be the solution of (36) defined on an interval [0, i e ndsurge], with xq(0) = x s i ng (X max ) and 
^o(^endsurge) = ^sing(7)- Now, to estimate the contraction, we linearize the layer system 

c(X-f(x)) 



Note that: 



d_ 
dx 



x 

X = 0. 

c(X-f(x)) 
f'{x) 



fix) 



(37a) 
(37b) 



/'(*) 

'f(xy 



\X=f(x) 

Hence, the first order coefficient of the contraction rate is estimated by: 



C 3 



iendsurge f'(x (t)) 



dt. 



f'(Mt)) 

Changing the variables and using (36) to express dt as dx/x' we get 

rW7) (f(x)) 2 g'(f(x))dx 



C3 



5 (x max ) f'(x)(f(x) + b ig (f(x)) + b 2 ) 



(38) 
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4.2 Canard phenomenon during the passage from pulsatility to surge 



While X increases from —7 to X max the point (x, y) may travel down along the left branch y = f(x) 
towards the fold (—Xf,—yf). Subsequently it either turns back and travels towards x S m g (X max ) 
along the left branch y = f(x) or jumps over to the other branch of y = f{x) to complete another 
pulse. The transition between these two possibilities is a canard phenomenon, corresponding to 
the travel along the middle branch of X. We will not describe this canard phenomenon in detail, 
restricting our attention to the computation of the maximal expansion. 

Proceeding as in §4.1 we consider the slow flow of the system (9), given by 

a x + aif(x) + a 2 + cX 
X ~ f(x~) ' (39a) 

X = -Y + g(X). (39b) 

and its solution (xo(t), Xo(t)) with initial conditions x(0) = —Xf and X(0) = Xf. The maximal 
amount of expansion occurs for trajectories traveling along the middle branch of y = f(x) from the 
vicinity of (—Xf, —yf) to (xf,yf) and is estimated by 

f ^max 

f'(x (t))dt 



where t max is the time when xo(t) reaches the upper fold (xf, yf), i.e. i max is defined by xo(t r 
Xf. Changing the variables, we get 



f(x (t))dt 



*' (f'(x)) 2 dx 

Xf a x + aif(x) + a 2 + cX (t(x)) 



where t(x) is defined by xo(t(x)) = x, with — Xf < x < Xf. Now sinc6 — <C X (t^ <C X m \ n wg 
have 

*f {f'{x)) 2 dx < f*f {f(x)fdx 

_ x a x + aif(x) + a 2 + cX (t(x)) J_ Xf ciqx + a\f(x) + a 2 - cXf ' 



Let 

j. f a x + aij{x) + a 2 - c^f 

It follows that the amount of contraction incurred due to the canard phenomenon is approximately 
equal to e Ci ^ e . 



r*t (f(x)fdx 

J- x , a x + aif(x) + a 2 - cX f ' 



4.3 Putting the pieces together 

In this section we put together all the phases of the dynamics, the surge, the pulsatility, the small 
oscillations, and the intermediate phases to get a transition around the entire cycle. We assume 
the 5 = 0(e), or more specifically, S < e. We introduce four sections of the flow corresponding to 
the different phases of the dynamics. 

Let rj be a small constant. We can now define the sections of the flow: 

Z™ = {(x,y,X,Y) : y = f(x f ) - r,}, 

S / = {(x,y,X,Y) : x = x f }, 

Vjsurge = y }XtY ) :X = X sing (X max ) + T)}, 

S endsu rg e = {{x ^ X ,Y) : X = X sing ( 7 ) " 7,}. 
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Figure 11: The four sections of the flow of system (1) represented in the (x,y)-plane. 

The transition from j] endsurge to S m involves a passage through a folded node point and hence 
there are canard type phenomena occurring there. However, the passage of the Regulator through 
this point is fast (0(1/8)). Hence we can assume that the trajectories stay away from the canards 
so that no significant expansion is present. 

If we restrict our attention to the interior of the sectors, staying away from the canards, then 
the transition from S m to £ surge involves no significant expansion. Close passage to canards can be 
due to the folded node near the upper fold or due to the canard phenomenon during the passage 
from pulsatility to surge described in §4.2. 

We can now quickly describe what happens if the trajectories stay away from the canards. 
Consider the region in S m consisting of the rotation sectors described in Theorem 2. The image 
of this region in E surge is contained in a compact subset of 2 surge . The image of this set by the 
transition from £ surge to vj endsur g e j s contained in a ball of radius 0(e~ C3 ^ s ), which we denote by B$. 
By adjusting a parameter, for example c, we can arrange that the image of B$ in S m is contained 
in the union of the rotation sectors. We first consider the case when the image of B$ in S in is not 
close to a canard. We can now consider the return map from j] endsurge to itself restricted to B$. 
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This transformation is well denned (maps Bg into itself) and is an exponential contraction. This 
proves the existence of a unique stable periodic orbit. 

We now consider the case when trajectories starting in £ surge pass near canards. One way that 
a canard segment can be involved in the recurrent dynamics is if the image of Bg in S m intersects 
two sectors (with e and 5 small enough, it cannot intersect more due to the estimate on the size 
of the sectors given in Theorem 2). There are now a few possibilities for how the trajectories can 
continue (see Figure 12). 

1. Trajectories that are not close to a canard transition to the pulsatility stage in the same as 
in the case when no passage near canards was involved. Near such trajectories no significant 
expansion is incurred. 

2. Trajectories that are close to canards, after passing through K2, return to S in through a 
segment of trajectory which looks like a canard cycle (with or without head). The maximal 
amount of expansion near such trajectories is incurred for maximal-like canards. The amount 
of expansion is estimated below. 

3. Small canards, that return to S-^ and then to EJf, without passing through S m . The amount 
of contraction near such trajectories is negligible. 

In both cases 2 and 3, X — Xf is already positive as the trajectory reaches X^"* so that a simple 
passage to pulsatility takes place following the return of the trajectory to £ m . 
















- ^ \ \ v^n 






2 / 
/ {' 


N 




y 


/ / \ 
/ ' >^ 
/ ' 


\ 

\ 

\ 

^- \ 




-1 


/ / 
/ / 
/ / 
/ / 
/ / 
// 


\ 

\ 

\ 

TT 

\ ' 

\ 


yin 


-2 


A 

1 \ 
1 \ 




\\ 

\ 

\ 

\ 

\ 

\ 


1 






\ 



3 1 — I , ^ , J — I 

-2 -1 1X2 



Figure 12: A "regular" canard and a small canard. 

The maximal expansion by this transition is e*^ 4 ^, where C4 is introduced in (40). To understand 
this estimate consider two points in S in which are endpoints of two trajectories starting in S m at 
two points very close to each other but on the opposite side of the maximal canard. The flow 
backwards in time can contract the distance by the maximal contraction, backward in time, along 
the middle part of the fast nullcline, which is bounded below by e~ C4 ^ e , with the constant C4 
computed analogously as C3, see §4.1. 

Another way that trajectories starting in £ surge may pass near a canard comes about by means of 
addition of a pulse at the end of pulsatility, or, in other words, by means of the canard phenomenon 
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described in §4.2. As described in §4.2 the maximal expansion is also bounded by e CA / £ . The 
cumulative effect of the expansion coming from the two sources is e 2C<4//e . 
Consequently, as long as 

^ < % (41) 

£ 

the return map from £ endsurge to itself is an exponential contraction of 0(e~ C3 ^ s+2C4 ^ £ ). Note that, 
for a fixed value of e, this condition is fulfilled for 5 small enough. In the estimate, we include the 
expansion incurred by the passage near both of the canard phenomena. Hence there exists a unique 
periodic orbit for every parameter value. This periodic orbit may contain a canard segment, which 
corresponds to a transition consisting of subtracting (respectively adding) a small oscillation and 
adding (respectively subtracting) a pulse at the beginning of the pulsatility stage, or it may include 
a canard segment corresponding to the addition of a pulse at the end of the pulsatility phase (the 
canard phenomenon described in §4.2). 



5 Numerical study 

We carried out numerical computation and continuation of periodic orbits in order to better describe 
the dynamics of our original four-dimensional system and visualize in more details the different 
transitions that shape the periodic orbits investigated in this work. We computed families of 
periodic orbits solutions of system (1) displaying pulses and surge, depending on various system 
parameters, using numerical continuation. We could then detect various canard-induced transitions 
affecting the number of pulses and the presence of a pause after the surge (see Figures 13, 14 and 15 
below). We also computed attracting and repelling slow manifolds, as well as secondary canards, 
near the folded node of system (18a)-(18c), which approximates the behavior of the full system 
during the pause and explain its small oscillations (see Figure 16 below). 

Systems with different time scales are well known to pose numerical problems because of their 
intrinsic stiffness, in particular when computing periodic orbits or, more generally, orbit segments 
[14, 15]. The use of numerical continuation in the context of slow-fast dynamical systems has 
significantly increased over the last two decades. First, in the classic framework of limit cycle 
continuation, where very sharp transitions such as canard explosions could be finely rendered. 
Second, and more recently [8], in the context of manifold computation, where slow manifolds were 
approximated by families of orbit segments computed by continuing a parametrized family of two- 
point boundary value problems (BVP). The combination of orthogonal collocation to compute 
an orbit segment solution of a BVP, with a predictor-corrector algorithm to move one or several 
parametrized conditions of the BVP, proved very efficient compared to shooting methods. Indeed, 
solutions of singularly perturbed ODEs display very sensitive dependence on initial conditions and 
parameter variations. In this context, orthogonal collocation gives a better approximation of such 
an orbit segment by distributing the error along the orbit instead of accumulating it at one end 
point as with shooting. Furthermore, the continuation algorithm gives a good rendering of the 
piece of manifold of interest, where orbit segments are distributed according to arc length, hence, 
accounting for the changes of local curvature of the manifold. In this way, one can integrate 
slow-fast ODEs with suitable boundary conditions using the BVP solver embedded in numerical 
continuation packages such as Auto [10]. We now illustrate the use of numerical continuation tools 
to investigate system (1). 
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5.1 Continuation of periodic orbits in parameter 02 

We start by a periodic orbit continuation that illustrates the various transitions, upon changes of 
parameter 02, that take place in between different parts of the typical periodic orbit of system (1) 
as shown in Figure 1. The other parameters are fixed at values previously fixed, that is, c = 0.69, 
ao = 1, a% = 0.02, d2 = 0.8, bi = 0, 62 = —0.8, A3 = —1, Ai = 1.5, ^3 = —1, fii = 4. The solution 
branch of periodic orbits is presented in Figure 13 where we choose to display on the vertical axis the 
maximum in y for each orbit as a measure of the solutions along the branch. The branch appears 
to be quite complicated with several rapid transitions that manifest themselves by quasi-vertical 
segments along the branch and that all have to do with canard trajectories. We identified two 
different types of transitions, affecting the periodic orbits at two different stages; before the surge, 
corresponding to the creation or annihilation of a pulse, and after the surge during the pause, 
corresponding to the transition of a small oscillation to a pulse. These transitions correspond 
to the canard phenomenon discussed in §4.2 and the canard phenomenon related to the small 
oscillations discussed in §3.1. Each transition takes place within an exponential small variation 
of the parameter, therefore thus corresponding to a quasi-vertical segment on the branch. We 
will now describe our numerical results, which reveal an intricate sequence of canard explosions 
corresponding to the two types of transitions. 

10.8 

max(i/) 
10.5 

10.2 



9.9 



0.75 0.77 0.79 0,2 0.81 

Figure 13: Family of periodic orbits solution of the original system (1) when 02 is varied. The 
vertical axis shows the maximum in y for each computed limit cycle along the branch. Eight orbits 
have been highlighted with black dots on the branch and given numbers from 1 to 8; they are shown 
in the two subsequent figures. 

To be more specific, we introduce a labeling scheme for periodic orbits of (1). We say that and 
orbit is of type (p, s) if it involves p pulses and s small oscillations. Our first transition can be 
described as (p, s) — > (p+ 1, s) (or (p, s) — )• (p— 1, s)), while the second one as (p, s) — > (p+ 1, s — 1) 
((p, s) —> (p — 1, s + 1)). We find two different scenarios in which both transitions occur. In the 
first scenario, they happen one after the other, which corresponds to two exponentially small bands 
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of parameter values with associated quasi-vertical segment of the branch, separated by an order 1 
interval of parameter values where the branch is, in comparison, quite flat. In the second scenario, 
they happen within the same exponentially small parameter variation. We isolate two sets of four 
orbits each along the branch, numbered 1 to 4 for the first set and 5 to 8 for the second, that 
undergo the first and the second transition, respectively. We now focus on each transition with the 
associated set of four chosen orbits. 

An example of the first scenario with both transitions is presented in Figure 14, where we show 
the time profile of y for the orbits 1 to 4 on the branch in Figure 13. This transition affects the 
number of small oscillations of the pause, this is why we enlarge each panel in the region of the 
pause and show the zoom in an inset; each panel is labeled with the number of the corresponding 
orbit in the solution branch. From orbit 1 to orbit 4, the pause gains one small oscillation, which 
corresponds to the loss of one pulse after the pause. However, one notices that there is very little 
difference between orbit 3 and orbit 4 in the pause. This is because this part corresponds to the 
second transition, where the change takes place before the surge with the appearance of one more 
pulse. In other words there are two canard explosions, which, using our labeling scheme, can be 
described as (8,3) — > (7,4) and (7,4) — > (8,4), giving a net result of an (8,3) — > (8,4) transition. 
Both transitions are canard- mediated, which one would see when plotting the orbits in the (x, y)- 
plane (see right panels of Figure 14). However, within this first scenario they are separated by an 
0(1) parameter interval. 

The second scenario is illustrated in Figure 15 where orbits 5 to 8 from the solution branch 
are represented in the phase plane (x,y). Here, both transitions seem to occur within the same 
exponentially small a2-variation. Hence, one has two canard explosions in this plane : the first 
one corresponds to the transformation of a pulse into an additional small oscillation on the pause 
(visible essentially from orbit 5 to 6), the second one corresponds to the gain of one pulse before the 
surge (visible essentially from orbit 6 to 8). Hence, starting from 7 pulses and 4 small oscillations 
in case 5, one obtains 7 pulses and 4 small oscillations in case 8. The whole transition can be 
described as (7,4) — > (6,5) — > (7,5). Note that the transitions in both scenarios are the same, 
the only difference is in the length of the parameter interval separating them. The two types of 
transitions are well explained by our theory, however we have no explanation for the very intriguing 
fact that, in the second scenario, there are two canard explosions occurring simultaneously. 

5.2 Computation of slow manifolds and secondary canards on the pause 

We now illustrate the change of small oscillations on the pause, upon parameter variation, by 
computing slow manifolds and secondary canards of system (18a)-(18c), which represents a good 
approximation of the full system (1) in the region of the pause. This system is three-dimensional 
and possesses a folded node for the parameter values we consider (see section 3.3). We computed 
slow manifolds and secondary canards near this folded node using the BVP strategy developed 
in [8, 9]. That is, we approximate the manifolds by a one-parameter family of orbit segments with 
initial conditions moving on a curve traced on the critical manifold C°, away from the fold J-, and 
end conditions restricted to a planar cross-section near the folded node. For simplicity, we take 
the one-dimensional manifold of initial conditions to be of the form C° n {x = xo} (where xq is 
chosen so that this line is at a large enough distance from the fold curve), and the two-dimensional 
manifold in which the end conditions lie to be of the form S enc [ = {X = X cn( j}, with Xo close to 
(So corresponds to a cross-section containing the folded node). The dimensions of these manifolds 
of boundary conditions are chosen so that the resulting BVP is well-posed. 
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Figure 14: First transition on the branch shown in Figure 13 upon variation of ai- From 1 to 4, 
the pause of the periodic attractor loses one small oscillations. We show the time profile of y to 
illustrate this transition on the left panels and the projection of the orbit onto the (a;, y)-plane on 
the right panels. 
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Figure 15: Second transition on the branch shown in Figure 13 upon variation of a<i- From 5 to 
8, the start of the surge undergoes a canard explosion. We show the time profile of y to illustrate 
this transition on the left panels and the projection of the orbit onto the (x, y)-plane on the right 
panels. 
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Figure 16: Attracting [Sf) and repelling {S r e ) slow manifolds of system (18) near its folded node. 
Panel (a) shows a representation of these manifolds together with three secondary canards in the 
3D phase space, in between the cross sections S_o.o25 := {X = —0.025} and So. 025 '■= {X = 0.025}. 
Panel (b) shows the intersection curves of the slow manifolds S£ and in £q := {X = 0}. 



In Figure 16, we show the result of these manifold computations. Panel (a) shows an attracting 
slow manifold (red) S° and a repelling one (blue) S£, computed in between sections £-0.025 and 
So. 025, together with three secondary canards (black curves) that correspond to transversal inter- 
sections between S® and S^; we also show the intersection curves (red and blue curves) of the slow 
manifolds with both cross sections. The spiralling behavior of the slow manifolds is typical of the 
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folded node scenario [27, 3]. In panel (b), we present the intersection curves of the slow manifolds 
with the cross section So that contains the folded node; once more, the figure is, as expected, very 
similar with previously computed slow manifolds in similar dynamical contexts [14, 27, 3, 8]. 

6 Conclusion 

In this paper we have studied the existence of MMOs in a system known as a phantom burster, 
consisting of two dimensional, unidirectionally coupled slow-fast oscillators or, in other words, one 
oscillator forcing the other by controlling the position of its null-cline. An additional feature of this 
system is the presence of three time scales; the dynamics of the forcing oscillator is slower than 
the dynamics of the one being forced. The orbits display alternatively three modes of dynamics: 
small oscillations near a fold, relaxation-type oscillations and a quasi steady state. Consequently, 
we needed to deal with quite complicated solutions, on the other hand the abundant structure of 
the system allowed us to obtain extensive results. 

We have obtained two main results, one of global nature, namely the existence of a unique 
attracting periodic solution, and one of local nature, analyzing secondary canards of a folded node 
with an additional slow time scale. The local result relies strongly on the additional slow time scale 
and gives an elegant and rigorous proof of the existence of secondary canards and of the sectors 
of rotation. The global result relies on the local result and on the existence of strong contraction 
during the quasi steady-state phase of the dynamics. This result is rather elementary in nature, 
but it yields a surprising conclusion: in certain regions of the parameter space, the transition from 
an MMO with n small oscillations to an MMO with n + 1 small oscillations is free of complicated 
dynamics, a unique stable periodic orbit exists through the canard transition. This way we have 
obtained a complete characterization of the dynamics for all values of the control parameter. 

We were able to obtain rather strong results due to the simplicity of the system in question, 
but we hope that some of the ideas and techniques can be extended to other situations, involving 
three time scale dynamics and/or having the phantom burster structure. One simple but powerful 
idea we used was to identify the phase of the dynamics of the slowest system (fast, slow, quasi 
steady-state) to design a reduced system fitting the part of the dynamics in question. This way 
we derived important reductions that greatly simplified the analysis. Clearly this technique will 
be applicable in other contexts of multiple time scale systems, even if there is two way coupling 
between the systems with different time scales. 

The continuation results of Section 5 gave a very nice and accurate illustration of our theo- 
retical results, as well as pointed to an intriguing phenomenon that we did not expect, namely 
a prediction of a double canard explosion. In this context it would be interesting to extend the 
continuation results to lower values of 8, which is a numerical challenge. A more detailed numerical 
and theoretical study of the double canard explosion will be a subject of future work. 

We would like to point out that the system we have studied has been used to model different 
modes of GnRH (Gonadotropic Releasing Hormone) secretion and transitions between them. The 
biological mechanisms underlying the transitions between the surge mode and the pulsatility mode 
in the physiological GnRH secretion pattern are still poorly understood. Therefore, our study may 
contribute to the development of tools and insights that can be used in the study of this very 
important problem. 
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